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ABSTRACT 

We present a 3-dimensional model of supernova remnants (SNRs) where the hydro dynamical evo- 
lution of the remnant is modeled consistently with nonlinear diffusive shock acceleration occuring at 
the outer blast wave. The model includes particle escape and diffusion outside of the forward shock, 
and particle interactions with arbitrary distributions of external ambient material, such as molecular 
clouds. We include synchrotron emission and cooling, brcmsstrahlung radiation, neutral pion produc- 
tion, inverse-Compton (IC), and Coulomb energy-loss. Boardband spectra have been calculated for 
typical parameters including dense regions of gas external to a 1000 year old SNR. In this paper, we 
describe the details of our model but do not attempt a detailed fit to any specific remnant. We also 
do not include magnetic field amplification (MEA) , even though this effect may be important in some 
young remnants. In this first presentation of the model we don't attempt a detailed fit to any specific 
remnant. Our aim is to develop a flexible platform, which can be generalized to include effects such as 
MFA, and which can be easily adapted to various SNR environments, including Type la SNRs, which 
explode in a constant density medium, and Type II SNRs, which explode in a pre-supcrnova wind. 
When applied to a specific SNR, our model will predict cosmic-ray spectra and multi-wavelength mor- 
phology in projected images for instruments with varying spatial and spectral resolutions. We show 
examples of these spectra and images and emphasize the importance of measurements in the hard X- 
ray, GeV, and TeV gamma-ray bands for investigating key ingredients in the acceleration mechanism, 
and for deducing whether or not TeV emission is produced by IC from electrons or pion-decay from 
protons. 

Subject headings: acceleration of particles — supernova remnants — cosmic rays — X-rays: general, 
gamma-ray 



1. INTRODUCTION 

Supernovae (SNe) are the only known sources capa- 
ble of providing the energy needed to power the bulk 
of the galactic cosmic rays (CRs) with energies below 
the s pectral feature called the "knee" around 3 x 10 15 eV 
(e.g.. lDrurvl[T983h . If SNe arc the main sources of Galac- 
tic CRs, the acceleration mechanism must be efficient 
so that > 10% of the total SN explo sion energy in our 
Galaxy ends up in cosmic rays (e.g., Hillas 2005). Ob- 
servational evidence that the outer blast wave shock ac- 
celerates electrons to ultra-relativistic e nergies in some 
young SNRs (e.g.. iKovama et all Il995h . and the exis- 
tence of a well-developed model of particle acceleration 
at shocks, i . e., diffusive shock acceleration (PSA) (e.g., 
D7urvlll983t iBlandford fc Eichleri fl987t Uones fc Ellison! 
19911 ) support the above contention. 

When confronting observations with theoretical 
models, however, there remain a number of important 
ambiguities and uncertainties from both the observa- 
tional and theoretical perspectives. Resolution of these 
ambiguities and uncertainties by new telescopes will be 
essential to claim evidence for the pion-decay feature in 
the GeV- TeV emission from SNRs. The Gamma-ray 
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Large Area Space Telescope (GLAST), to be launched 
in 2008, will probe this crucial energy range with un- 
precedented sensitivity and resolution. 

Fundamental questions for CR origin also concern 
the spectral shape and maximum ion energy a given 
SNR can produce. Electron energy spectra inferred from 
young SNRs vary and can be substantially harder than 
CR electron spectra observed at Earth, even after correc- 
tion f or propagation in the galaxy (e.g.. lBerezhko fc V51kl 
12001 . The maximum CR ion energy SNRs actually pro- 
duce will remain uncertain until a firm identification of 
pion-decay emission is obtained and gamma-ray emission 
is detected past a few 100 TeV, the maximum possible 
electron energy in SNRs. 

There remain other basic questions concerning the 
DSA mechanism. For instance, is DSA efficient enough 
for nonlinear effects, such as shock smoothing and mag- 
netic field amplification, to become important in young 
SNRs? How does particle injection occur and how does 
injection and acceleration vary between electrons and 
protons? While the galactic CR electron-to-proton ratio, 
(e/p) re i, of 0.01-0.0025 observed at Earth at relativistic 
energies is often used to constrain the ratio in SNRs, this 
ratio has not been observed outside of the heliosphere. 1 

1 We note that while energetic electrons and protons are observed 
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The (e/p)rei ratio is crucial in deciding whether the 7- 
ray emission from different SNRs, or observed in different 
parts of an individual SNR, is of hadronic or leptonic ori- 
gin. 

The recent discovery of spatially t hin, hard X- 
ray f i laments in some youn g SNRs (e.g., iBamba et alJ 
120031 ; lUchiyama et al.l l2007f) supports previous sugges- 



tions (e.g.. iCowsik fc Sarkail fl980l : iBell fc Lucekl 12001 



iRevnolds fc Ellison! 1 19921 ) that the particle acceleration 
process can amplify the ambient magnetic field by large 
factors. If magnetic field amplifica tion in PSA is as larg e 
as now appears to be the case (e.g.. lBerezhko et aI1l2003l ). 
it will have far-reaching consequences not only for under- 
standing the origin of Galactic CRs, but for interpreting 
synchrotron emission from shocks throughout the uni- 
verse. Since shocks and related superthermal particle 
populations exist in diverse environments, the knowledge 
gained from studying SNRs will have wide applicability. 

The advent of new space- and ground-based tele- 
scopes will result in observations of SNRs at many 
different wavelengths with greatly improved sensitivity 
and resolution. It is even conceivable that features in 
the CR spectrum observed at Earth might be associ- 
ated with nearby SNRs with future observations (e.g., 
lErlvkin fc Wolfendaldfl999l : iKobavashi et~aHl200ll) . 

In order to take full advantage of current and fu- 
ture observations, and to improve our understanding of 
the DSA mechanism, the data must be analyzed with 
consistent, broadband photon emission models including 
nonlinear effects. This has prompted us to develop a 
three-dimensional model of young SNRs where the evo- 
lution of the remnant is coupled to nonlinear diffusive 
shock acceleration (NL-DSA) (e.g., lEllison et all 120041 ; 
lEllison fc Cassam-Cheiial2005| ). in an environment with 
an arbitrary mass distribution. We focus on radiation 
from CR electrons and protons and leave the modeling 
of heavier ions for future work. In this preliminary study, 
we also ignore other possible acceleration processes, most 
notably second-order stochastic acceleration, and do not 
include magnetic field amplification. 

We believe our work is a significant advance over pre- 
vious work for several reasons. Of particular importance 
is that we include "escaping" particles self-consistcntly. 
In NL-DSA, a sizeable fraction of the SN explosion en- 
ergy can be put into very energetic CRs that escape 
the forward shock and stream into the surrounding ISM. 
These particles will produce detectable radiation if they 
interact with dense, external material. Another advan- 
tage is that we have a "coherent" model, easily expand- 
able to include more complex effects, where the var- 
ious environmental and theoretical parameters can be 
straightforwardly varied and the resulting radiation can 
be compared directly with observations. This is impor- 
tant because all SNe and SNRs are different and com- 
plex with many poorly constrained parameters. It is 
essential that the underlying theory consistently model 
broad-band emission from radio to TeV 7-rays taking 
into account individual characteristics of the remnants 
and their environments. 

from solar flares and at low Mach number heliosphcric shocks, these 
observations provide limited help for understanding the high Mach 
number shocks expected in young SNRs and other astrophysical 
sourses where a large fraction of the shock energy is put into rela- 
tivistic particles. 



In Sections [2] and [3] we give a brief general descrip- 
tion of nonlinear diffusive shock acceleration and describe 
the environmental and model parameters required for a 
hydrodynamical solution. We place a time-dependent, 
spherically symmetric, hydrodynamic calculation of a 
SNR, including NL-DSA, in a three-dimensional box con- 
sisting of 51 x 51 x 51 cells. 2 The energetic particles pro- 
duced by the outer blast wave shock propagate through 
the simulation box where they interact with an arbitrary 
distribution of matter placed external to the outer shock. 
The energetic particles in the box, including those within 
the SNR, suffer energy losses and produce broad-band 
continuum emission spectra by interacting with the mag- 
netic field, photon field, and matter density of each cell. 
In Section [4] we show some examples including line-of- 
sight projections of the emitted radiation which are suit- 
able for comparison with observations. 

There are a number of young SNRs under 
active i nvestigation, in clu ding SN R RX J1713 
' e.g.. lAharonian et al.l 120061; lUchiyama et all 
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Rho et all 12002ft. IC 443 
VERITAS Collaboration: T B. Humenskvl 
and W 28 (|Aharonian et alJ 12008ft " However, here we 
concentrate on a general study using various parameters 
typical of young, shell Type la SNRs and leave detailed 
modeling of individual remnants for future work. 

2. DIFFUSIVE SHOCK ACCELERATION IN SNRS 
2.1. The Diffusive Shock Acceleration Theory 

In the test-particle approximation, diffusive shock 
acceleration produces superthermal particles with a 
power law distribution where the power-law index de- 
pends only on the shock compression ratio, i.e., f(p) oc 
where a = 3rTp/(rTp — 1), r^p is the test-particle 
shock compression ratio, p is the particle momentum, 
and f (p) is the phase-space dist ribution function (see 
lDrurvlll983t iBlandford fc Eichlerlll987l and references 
therein). This test-particle result holds as long as the 
pressure exerted by the accelerated particles (i.e., cosmic 
rays), P cr , is small compared to the far upstream mo- 
mentum flux, pqUq (pq is the unshocked density and uq 
is the unmodified shock speed). There is considerable 
observational evidence, however, that DSA is intrinsi- 
cally efficient and shocks with high sonic Mach numbers 
Ms > 10 are expected to accelerate particles efficiently 
enough that P cr ~ p Q u^. In this case, the pressure in ac- 
celerated particles feeds back on the shock structure in a 
strongly nonlinear fashion (e.g.. iJones fc Ellisonl[l99l . 

In NL-DSA, the following effects become important: 
(i) a precursor is formed upstream of the viscous sub- 
shock with a length scale comparable to the diffusion 
length of the highest momentum particles the shock 
produces. In the shock reference frame, the incoming 
plasma is decelerated and heated in the precursor before 
it reaches the subshock; (ii) the production of relativistic 
particles, and the escape of some fraction of the highest 
energy particles from the precursor, soften the equation 
of state of the plasma, making the plasma more com- 
pressible and allowing the overall shock compression ra- 

2 The resolution of the 3-D box is, of course, adjustable and 
limited only by computational considerations. 
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tio to increase, i.e., r to t > ^tp! (iii) the simple power 
law of the test-particle approximation is replaced by a 
concave spectrum at superthermal energies. The spec- 
trum is softer than the test-particle power law for low 
momentum particles and harder for high momentum par- 
ticles; and (iv) the weak subshock has a compression ra- 
tio r su b < ttp so that the shocked plasma has a lower 
temperature than would be the case in the t est-particle 
approximation (see iBerezhko fc Ellis"onlll999l and refer- 
ences therein for detailed discussions of these effects). 

The modification of the equation-of-state by the pro- 
duction of relativistic particles and the escaping energy 
flux in NL-DSA, influences the evolution of the SNR and 
numerical approaches have been developed to describe 
this process ('e.g.. |Berezhko. Elshin. fe Ksenofontovl 
119961 Ellison. Decourchelle. fc Ballet! \206W . Here we 
generalize the basic NL-SNR model by including CR 
propagation within the remnant and, most importantly, 
in surrounding material using a three-dimensional simu- 
lation. The escaping particle flux is expected to dominate 
interactions outside of the SNR blast wave. 

2.2. CR-Hydro Simulation 

We calculate the hydrodynamic evolution of 
a SNR with a sp herically symmetric model de- 
scribed in detail in I Ellison et al.l (|2007f ) and refer- 
ences therein (see Fig. [TJ). The model couples ef- 
ficient DSA to t he hydrodynamics using the semi- 
anal ytic model of iBlasi. Gabici. fc Vannonil (|2005f ) (see 
also lAmato fc Blas i 2005, 2 0061) . Given an i njection pa- 
rame ter, Xinj (this is £ in equation (25) in IBlasi et all 
I2005T ). the semi-analytic model calculates the full proton 
distribution function f p (p) at each time-step of the hydro 
simulation, along with the overall shock compression ra- 
tio, rtotj an d the subshock compression ratio, r su b- The 
hydro provides the required input for the semi-analytic 
calculation, i.e., the shock speed, shock radius, ambi- 
ent density and temperature, and the ambient magnetic 
field, and f p (p) reflects the nonlinear effects from effi- 
cient acceleration. The coupling between the hydro and 
NL-DSA is accomplished by using f p (p), and the escap- 
ing particle flux, to calculate an effective ratio of specific 
heats which is then used in the hydrodynamic equations. 
The electron spectrum, f e (p), is determined from f p (p) 
with two additional parameters, the electron-to-proton 
ratio at relativistic energies, (e/p) Te i, and the tempera- 
ture ratio immedi ately behind the shock, (T e /T p ) (see 
lEllison et al.ll2004 for a full discussion). 

In this paper, we only consider Type la supernovae 
with no pre-SN wind. We also ignore any CR production 
that might occur at the reverse shock. Both of these re- 
strictions are for clarity and our model can be applied to 
Type II SNe with winds and can calculate particle heat- 
ing and acceleration at reverse shocks. The parameters 
controlling our results fall into two catagories: Environ- 
mental parameters and model parameters. These arc 
listed in the following sections with either default values 
or the range of values used for our examples. 

2.2.1. Environment Parameters 

The environmental parameters include: (i) the SN 
explosion energy, E$n = 10 51 erg, (ii) the ejecta mass, 
M e j = 1.4 M Q , (iii) the distance to the SNR, £> SNR = 
lkpc, (iv) the age of the SNR, £ S nr = 1000 yr, (v) the 



ISM proton number density, n p — 0.1, 1, or 10 cm" 3 , 
(vi) the proton number density in the molecular cloud 
if present nyic = 10 3 cm~ 3 , (vii) the ambient, i.e., un- 
shocked, magnetic field, £>ism = 3 fiG, and (viii) the am- 
bient proton temperature, T p = 10 4 K. The quantities 
n p , -BisM; and T p are assumed to be constant in the re- 
gion outside of the forward shock. 



2.2.2. Model Parameters 

The model parameters used in this simulation are: 
(i) an exponential ejecta density profile applicable to 
Type la SNe, (ii) the acceleration efficiency for DSA, 
face = -E-cr/^sn, where we consider two possibilities: 
the test-particle case where 1% of the total SN explo- 
sion energy is put into CR energy, Eqr, during the 
1000 yr evolution of the SNR, and nonlinear DSA, where 
75% of the SN explosion energy is put into CRs during 
1000 yr, 3 (iii) the electron to proton ratio at relativis- 
tic energies, (e/p) rc \ = 0.01, (iv) the electron to proton 
temperature ratio immediately behind the forward shock, 
(T e /T p ) = 1, (v) the cutoff index for the shape of particle 
spectra near -E ma x, ct cu t = 1) ( v t ne number of gyroradii 
in a mean free path, 77 m f p = 1, (vii) the fraction of the 
forward shock radius, / s k = 0.05, used to truncate DSA, 5 
(viii) the number of shells between the forward shock 
and the contact discontinuity at the end of the simula- 
tion, -/Vsheii = 20, and (ix) the diffusive time step interval, 
£stcp = 10 yr. All of these pa rameters, except bmc and 



t s tep, a re d escribed in detail in Ellison fc Cassam-Chena'3 
(|2005h and lEllison ei~aLl J20Q3). 

The geometry of the magnetic field that is used 
as input to the DSA calculation and to calculate the 
synchrotron emission is not described explicitly in the 
CR-hydro simulation. Instead, it is assumed that 
the field immediately upstre am from the FS , Bq = 
-Bism, is turbulent and, as in iVolk et all j2002), we set 
the immediate downstream compressed field to B2 = 
B ^/l/3 + 2rf ot /3, where r to t is the overall shock com- 
pression ratio. The magnitude of the shocked field 
ev olves as the density of the p l asma changes, as described 
in lEllison fc Cassam-Chenaa (|2005| ). and the magnetic 
pressure is included in the hydrodynamics, although it is 
insignificant for the results we show here. An important 
limitation of our current model is that we do not include 
self-generated magnetic turbulence or magnetic field am- 
plification. Magnetic field amplification i s only now being 
studi e d in nonlinear calculat ions (e.g., lAmato fc Blasil 
120061 : IVladimirov et al.1 120061 ) and we leave implemen- 
tation of this important aspect of DSA for future work. 
We also neglec t other wave-par ticle effects, such as wave- 
damping (e.g.. lPohl et aLlfeOOSf L and simply assume that 
the shocked field is turbulent enough for Bohm diffusion 
to occur with a background field that is compressed at 
the shock and evolves adiabatically behind the shock. 

3 These percentages include CRs that escape upstream from the 
forward shock during the SNR evolution. 

4 This parameter is discussed more fully in Section 13. 1 1 below. 

5 The maximum proton energy produced by the shock, -B max , 
is determined by either the finite shock age, tsNR, or the fi- 
nite size of the shock, whichever occurs first. Our choice of 
/sk = 0-05 is arbitrary but is c onsistent with previous work (e.g., 
Ellison & Cassam-Chenai 2005). For this particular f^, E max IS 
determined by the finite shock size in all of our examples. 
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Fig. 1. — Cross-section of the 3-D boxcl simulation (not to scale). 
For the results shown in this paper, the box is divided into 51 X 51 X 
51 cells. Only a few representative SNR shells are shown in this 
sketch. The actual number of spherically symmetric shells between 
the forward shock (FS) and contact discontinuity(CD) increases 
with time and equals 20 at the end of the simulation. We show the 
'molecular cloud' discuss ed in Section 14.41 The 'molecular cloud' 
shell discussed in Section |4. 3. 21 is not shown for clarity. 

2.3. Model Geometry and Simulation Method 

We treat the SNR hydrodynamics in 1-D by assum- 
ing a spherically symmetric structure for the region of 
the remnant between the forward shock (FS) and the 
contact discontinuity (CD). The ma in generalizat i on we 
have made to the CR-hydro model of lEllison et al.l (|2007f ) 
is to imbed the SNR in a fully 3-D astrophysical environ- 
ment where CRs accelerated by the remnant propagate 
and interact with ambient material. A cross-section of 
the 3-D simulation box is shown in Fig. [TJ Spatially de- 
pendent environmental aspects, like matter density in a 
molecular cloud, magnetic field strengths, and the mag- 
netic turbulence spectrum are all defined and stored in 
3-D simulation cells. 

The various interactions and photon emission pro- 
cesses are computed throughout the simulation box so 
that mult-wavelength spectra and projected morphology 
are obtained. 

The temporal sequence of the evolving SNR is as fol- 
lows: For a SNR of a given age, tsNR, we divide its life 
span into iV s heii epochs. During each epoch, the forward 
shock propagates into the ambient medium and new CRs 
(and shock-heated ISM plasma) are produced. A new 
spherical shell containing the shocked thermal plasma 
and new CRs is created. In every subsequent epoch, 
this spherical shell of material evolves while another new 
shell is produced. In this way, an "onion skin" struc- 
ture is formed of shells containing CRs of various ages 
(see Fig. [T]). The evolution of the shells includes the 
hydrodynamics (i.e., adiabatic effects), changes in the 
assumed frozen-in magnetic field, and losses from radia- 
tion and Coulomb processes for electrons. Spatial diffu- 
sion of CRs, magnetic field evolution in the shells, and 
fast synchrotron losses for electrons are treated using a 
finer timescale through further dividing each epoch into a 
number of time steps t st0 p = 10 yr. As the local magnetic 
field evolves in the shocked material, the local diffusion 



coefficient is modified accordingly. 

As mentioned in footnote El the maximum CR en- 
ergy in our examples is determined by the finite shock 
size. Particles that reach this energy escape and, for ef- 
ficient DSA, carry away a sizable fraction of the total 
energy flux. 

For each epoch, the CR-hydro simulation determines 
the escaping flux and maximum CR energy, E nmx , for 
electrons and protons in the outermost shell immedi- 
ately behind the FS where CR acceleration is taking 
place. These particles are added to the simulation box 
in a spherical shell immediately in front of the FS. While 
the precise energy distribution of the escaped particles 
is still largely unknown (the shape is not determined by 
the CR-hydro model), we assume the escaped CRs have 
a Gaussian distribution in momentum-space centered at 
-E-max and normalized to the total escaped flux (E max and 
the tot al escaped flux are determin ed by the CR-hydro 
code) dZirakashvili fc Ptuskin] I2008D . The width of the 
Gaussian is determined by fitting the high-energy spec- 
tral cut-off around E max of the newly accelerated CRs. 
The width of this cut-off depends on our model parame- 
ter, a cut . 

As time progresses, the energetic electrons and pro- 
tons diffuse in both the SNR shells and in the exter- 
nal material with momentum-dependent diffusion coeffi- 
cients described in the next section. As the CRs diffuse, 
they interact with the astrophysical environment, such as 
the cosmic microwave background (CMB) radiation or a 
molecular cloud, and the photon emissivity is recorded 
as a 3-D map for later analysis. 

3. DIFFUSION AND INTERACTION PROCESSES 
3.1. Diffusion 

In the simulation, particles spatially diffuse in two 
distinct regions; the volume inside the shocked SNR 
shells and the ambient ISM outside of the FS. For the 
shocked material, we assume Bohm diffusion while for 
the unshocked ISM we assume much weaker diffusion. 
For this study, we assume Kolmogorov turbulence domi- 
nates outside of the SNR but other forms could be used 
instead. Inside and outside of the SNR, we assume the 
turbulence is strong enough to ensure isotropic diffusion 
over length scales < 1 pc. 

If A is the scattering mean free path, the diffusion 
coefficient can be written as: 



D = -Xv , 

3 ' 



(1) 



or, if we assume A is proportional to some power of the 
gyroradius, 

D = D a [3(r a /r g0 ) s , (2) 

where v is the particle speed, (3 = v/c, r g = pc/{eB) is 
the gyroradius in cgs units, r g o is some constant reference 
length, Do = ?7mfp^goc/3 is a normalization constant, and 
s depends on the magnetic turbulence spectrum. 
For Bohm diffusion, r/ m f p = s = 1 and 



D, 



pc 



3 \eB(r,t) 



(3) 



Bohm diffusion is assumed throughout the shocked gas 
and the magnetic field in a particular shell, i?(r,i), de- 
pends on the location of the shell and its age. We assume 
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Fig. 2. — Momentum-dependence of spatial diffusion coefficients 
for CR protons. Bohm diffusion (dashed curve) is implemented for 
the space inside the SNR shells (-Db); The Kolmogorov spectrum 
(solid curve) is employed for the space outside the shells for CR 
diffusion in the ISM {D Kol ). A field B = 128 fiG is used for the 
Db plot. The normaliza tion of DkoI is taken from the calculation 
for the Galactic ridge by[Ptuskin ct al. (2 006T) . 

B(r,t) is 'fr ozen-in' and the details of t he fie ld evolution 
arc given m Ellison fc Cassam-Chen"al (|2OO50 . 

For the volume outside of the FS, Kolmogorov tur- 
bulence is assumed (s = 1/3) and the no rmalization of 
the di ffusion coefficient -DkoI 1S taken from lPtuskin et alj 
(|2006[ ). a value determined to reproduce the observed CR 
spectra at Earth, i.e., 



£>Koi = 0.25/3 



10GV 



1/3 kpc^ 
Myr 



(4) 



where R = pc/(Ze) is the magnetic rigidity. The dif- 
fusion coefficients are shown in Fig. [2] and, as expected, 
Db -C -DkoI since the self-generated turbulence in the 
shocked material is far stronger than turbulence in the 
relatively undisturbed ISM. 67 

Simple diffusion of CR particles is incorporated in 
the simulation in a discrctized manner. In each time step 
and each spatial grid in the 3-D simulation box, particles 
are exchanged between the adjacent boxels according to 
the particle momentum, location, and density gradient. 
The particle's location determines which diffusion coef- 
ficient is used, and the simulation resolution is mainly 
determined by the boxel size and time step, i s tcp, which 
are user-tunable. 

3.2. Interaction Processes 

The CR interaction processes considered include 
synchrotron radiation, bremsstrahlung, invcrse-Compton 
scattering, and neutral pion decay. Energy changes from 
adiabatic effects and radiation, as well as Coulomb en- 
ergy losses, are also included. All of these processes 

6 We make no attempt to self-consistently calculate the turbu- 
lence generated by CRs as they escape from the SNR and stream 
through the ISM. 

7 We note that Eq|4]is significantly differe nt from one assumed 
in a recent paper Gabici & Aharonian (2007) where a strategy to 
search for "PeV accelerators" in SNRs is discussed. The difference 
is due to their assumption that generation of plasma waves can 
suppress the diffusion coefficient by an order of magnitude relative 
to that for Galactic cosmic rays. 



arc treated in a fully space- and time-dependent fashion 
where the evolution of relevant parameters, such as the 
magnetic field and shell densities, are taken into account 
in each time step and boxcl. 

The details o f the radia tion processes can b e found 
in lSturner et all (|1997f ) and iBaring et all (|1999f ) but we 
note that, for IC emission, we only consider CR electrons 
colliding with a monoenergetic and isotropic photon field 
with an average energy density equal to that of the CMB 
field. 

For hadronic interactions we employ the latest 
parametric proton -proton (p-p) model developed by 
iKamae et al.1 ((2006). In this model, the total inclusive in- 
elastic p-p cross-section includes the non-diffractive (with 
Feymann scaling violation) and diffractivc components, 
plus the A(1232) and i?es(1600) resonance-excitation 
contributions important in the lOMeV to 1 GeV range. 
This model alone can account for ~ 20% of the GeV 7 
ray excess between the EGRET Galactic diffuse spec- 
trum and previous model prediction u sing proton data 
in the solar system (|Hunter et al.lfl997t h 

3.2.1. Coulomb Losses 

Coulomb losses for superthermal electrons are 
included in our m odel using equation (10) from 
ISturner et all ()1997f) . i.e., 
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A(t)rtSNR?7He 



m)-i>'(t)}, (5) 



where ngNR is the proton number density in a shocked 
shell, (3 e — v e /c is the electron /3, and t is the time. 
The d efinitions of the other terms are in ISturner et al.l 
(|1997f ) but are not important for our discussion here. 
Equation (JSJ) shows that Coulomb losses increase for large 
ambient densities and low electron speeds. As a shell 
of shocked material ages, Coulomb losses cause the low 
energy part of the superthermal electron distribution to 
become depleted, as indicated in Fig. [3J In all cases, 
Coulomb losses are insignifcant for protons. 

4. RESULTS 

In this initial presentation of our 3-D simulation, 
we show results for a set of generic Type la SNR mod- 
els where we vary the acceleration efficiency for DSA, 
£acc = Ecr/Esn, and the ambient proton number den- 
sity, n p (we assume the ISM is made of hydrogen). All of 
the other environmental and model parameters are kept 
constant with the values given in Sections l2.2.1l and l2.2.2l 

The values of e acc and rip for the four models we 
compute are given in Table [1] We have included the 
injection efficiency, Xinj, m Tablc[T] where Xinj determines 
Ecr/Esn and in practice we vary Xinj until we obtain the 
desired acceleration efficiency. We also show the fraction 
of -Esn that is in escaping particles, e OS c- Model A is 
used as a reference for the other three and, for all of the 
models, the duration of each epoch is 50 yr so we have 
20 shells in the SNR when it is 1000 yr old. 

4.1. Electron and Proton Spectra 

In Fig. [3] we show electron and proton phase-space 
distributions for the four models listed in Table [TJ We 
plot p 4 /(p) to emphasize the spectral curvature at rela- 
tivistic energies and the spectra are integrated over the 



6 Lee et al. 



TABLE 1 

SNR Model Parameters & Results 



Model Eacc Xinj rip i?FS a rtot b Pmax C e asc d -FkcV/ToV C ^kcV/TcV.MC 

(cm- 3 ) (pc) (10 4 m p c) 

A 0.75 3.70 1.0 4.67 10.66 1.37 0.27 1.24 0.97 

B 0.01 4.27 1.0 5.14 4.04 1.77 0.0015 74.0 68.5 

C 0.75 3.43 0.1 6.73 10.35 3.29 0.25 5.66 4.59 

D 0.75 3.77 10.0 3.06 11.89 0.53 0.28 0.19 0.16 



a Radius of the forward shock at the end of the simulation 

b Total compression ratio at the end of the simulation 

c The maximum momentum of CR protons at the end of the simulation 

d Fraction of i?SN carried away by the escaped protons at the end of the simulation. Note that e acc includes this fraction. 

c Energy flux ratio between emission at energies of 3 keV and 1 TeV at the end of the simulation. The fluxes are integrated over 

energy bands with widths of 1/10 of the central energies, and over the entire source volume. 

f Energy flux ratio as above, but now also including emission from the shell molecular cloud described in section 14.3.21 



entire shocked region between the FS and CD at the end 
of the simulation. 

In the top panel, we compare Models A (efficient 
NL-DSA; bold lines) and B (inefficient DSA or TP ac- 
celeration; thin lines). The test-particle model shows 
flat electron and proton spectra at relativistic energies 
[/(p) K P~ 4 ] with considerably lower fluxes at relativis- 
tic energies than Model A. The 'thermal' portions of the 
spectra show that the TP shock produces higher temper- 
atures than in Model A, a characteristic feature of NL- 
DSA. The structure seen in the 'thermal' portions of the 
spectra comes about because these spectra are summed 
over the various shells and the ones produced early on 
have less efficient DSA and have a higher temperature 
than later shells. 

In the middle and lower panels of Fig. [3] we keep 
£acc = 0.75 but vary n p ; n p = 0.1 cm -3 in the middle 
panel and n p = 10 cm -3 in the lower panel. The im- 
portant points for this comparison arc: (i) the CR flux 
at relativistic energies scales approximately as n p , as ex- 
pected, (ii) the maximum prot on momentum scal es in- 
versely as n p (see, for example. iBaring et al.lll999l ). (iii) 
the electron cutoff energy also scales inversely as n p but 
is influenced by radiation losses and the dependence is 
weaker than for protons, and (iv) the shocked tempera- 
ture scales inversely as n p , although this is not immedi- 
ately clear from the figures since the 'thermal' portions 
of the distributions are made up of contributions from a 
range of temperatures and densities. 

Of course, other aspects of the hydrodynamics de- 
pend strongly on n p . The radius of the SNR at £snr = 
1000 yr is considerably greater for n p = 0.1 cm -3 (i?FS = 
6.7pc) than for n p = 10cm -3 (Rfs = 3.1 pc). It is 
also expected that the FS will weaken faster with time 
for a denser upstream medium. However the strength, 
in terms of the efficiency of NL-DSA, also depends on 
the magnetic field and for the parameters used here, 
Model D has a larger compression ratio at isNR = 1000 yr 
(r tot = 11.9) than Model A (r tot = 10.7). 

Coulomb losses also increase as n p increases and the 
dip which appears just above the thermal peak in the 
Model D electron spectrum (light dashed curve in lower 
panel) reflects Coulomb losses experienced by the su- 
perthermal electrons as they collide with the shocked 
thermal gas. Coulomb losses can be expected to be more 
pronounced in NL models because the larger compression 
ratio results in a larger post-shock density. 



4.2. Spatial Variation and Escaping Flux 

At any given time, the spatial variation of the CR 
spectrum can be calculated. Fig. 2] shows CR spectra of 
Model A at three different locations: (i) just behind the 
forward shock (solid lines) ; (ii) mid- way between the for- 
ward shock and the contact discontinuity (dashed lines); 
and (iii) at a distance of d = 9 pc from the center of the 
SNR which is approximately 2 — 6pc beyond the FS, de- 
pending on the model (dotted lines). The heavy-weight 
curves are protons and the light-weight curves are elec- 
trons. 

Compared with the freshly accelerated electrons at 
location (i), many of the highest energy electrons are lost 
in the mid-point location (ii) mainly due to synchrotron 
losses with a small contribution from adiabatic losses. 
The protons show a smaller change which is due to adi- 
abatic losses only. 

At location (iii), only those CRs that escaped from 
the shock are present and their spectra lack a low energy 
component since low energy CRs remain trapped in the 
remnant. The hardness of the spectra at 9pc from the 
center of the SNR reflects the strong momentum depen- 
dence of the escape probability and the spatial diffusion 
coefficients. The escape probability from the SNR in- 
creases with energy and high-energy CRs diffuse faster 
in the ISM. 

4.3. Boardband Photon Spectrum 

Once the particle spectra are determined, the pho- 
ton emission can be calculated throughout the simulation 
box for arbitrary 3-dimensional distributions of matter 
and ambient photon fields. 

We consider two simple matter distributions (i.e., 
'molecular clouds') outside of the FS. The first is a spher- 
ical shell, concentric with the SNR where the inner and 
outer radii are equal to 9 and 10 pc, respectively. The 
second is a hemisphere centered at one side of the simula- 
tion box with radius = 3.2 pc (see Fig. [1]). In both cases, 
the proton number density in the 'molecular cloud' is 
Hmc = 10 cm" 3 and the magnetic field is -Bism = 3/xG, 
the same field as in the ISM. The entire simulation box 
is 20 pc on a side and is divided into 51 x 51 x 51 boxels. 
The density in the ISM between the molecular clouds 
and the FS is n p and the photon field throughout the 
simulation box is the uniform CMB field for all models. 

4.3.1. Emission from the SNR Shells 
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Fig. 3. — CR spectra integrated over the whole SNR FS-CD 
region at 1000 yr, plotted as p f(p), where f(p) is the phase space 
distribution function. In all panels, solid curves are protons and 
dashed curves are electrons. Also in all panels, the heavy curves are 
from Model A and these are compared to thin curves for Model B 
(top panel), Model C (middle panel), and Model D (bottom panel). 
Parameters for the various models are given in Table [T] 



Fig. [5] shows the broadband photon emission for 
Models A to D integrated over the shocked SNR shells 
between the FS and CD. The bottom panel shows the 
total spectra while the upper three panels show the in- 
dividual components from 7r°-decay (solid) , IC (dashed) , 
synchrotron (dash-dotted), and bremsstrahlung (dotted) 
compared with Model A. Emission from CRs outside of 
the SNR is not shown. 

In Model A, the photon flux in the radio to X-ray 
energy range is dominated by synchrotron emission up to 
~ 100 keV. The second largest contribution is from ther- 
mal bremsstrahlung which dominates between ~ 100 keV 
and ~ 50MeV. Between ~ 50MeV and ~ lOGcV, pion- 
decay and IC compete. Beyond ~ 10 GeV, the emission 
is dominated by IC. 
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Fig. 4. — Spatial variation of CR momentum distribution for 
Model A at t = lOOOyr. The distribution is plotted at three loca- 
tions: (i) just behind the FS (solid); (ii) mid-way between the FS 
and the CD (dashed); and (iii) at a distance of 9pc away from the 
SNR center (dash-dot). Thick lines are for protons and thin lines 
are for electrons. 



As seen in the A vs. B comparison panel, thermal 
bremsstrahlung plays an important role in the TP Model 
B and dominates synchrotron emission in the entire X- 
ray energy band. Thermal bremsstrahlung is strong in 
the TP model because the shocked temperatures are con- 
siderably higher than those in efficient DSA. The emis- 
sion from synchrotron, IC, and pion-decay are all weak 
in the TP case as expected. 

In the three NL Models A, C, and D, the acceleration 
efficiency is set at e acc = 0.75, but the ambient density 
is varied with n. p = 1, 0.1 and 10 cm -3 respectively. In 
the X-ray band, the thermal bremsstrahlung scales ap- 
proximately as 7ip and dominates synchrotron in Model 
D, where n p = 10cm -3 . 

Above ~ 100 MeV, the competition is mainly be- 
tween IC and pion-decay but bremsstrahlung is also im- 
portant for n p = 10 cm" 3 . For Model C (n p = 0.1cm -3 ), 
both pion-decay and bremsstrahlung are suppressed rel- 
ative to IC. For Model D, pion-decay dominates until 
near the maximum energies where bremsstrahlung be- 
comes comparable. 

4.3.2. Emission from a Shell 'Molecular Cloud' 

We first consider the shell of external material cen- 
tered with the SNR. Protons and electrons which have 
sufficiently high energy and, therefore, long diffusion 
lengths can escape from the FS and enter the ISM. These 
CRs also interact with the ambient ISM material of den- 
sity rip and the CMB radiation. 

In Fig.[6]we show the photon spectra from the molec- 
ular cloud shell for Models A-D in the same represen- 
tation as Fig. [5] but now integrated over the molecular 
cloudvohtme, all calculated at 4snr = 1000 yr. As ex- 
pected, these spectra are considerably harder than their 
counterparts inside the remnant. The escape of CRs 
from the forward shock during acceleration depends on 
the diffusion coefficient and the strong momentum de- 
pendence of the Bohm diffusion coefficient, Db, favors 
the escape of the highest energy particles. Once in the 
ISM, the rclativistic CRs diffuse with a diffusion coeffi- 
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Fig. 5. — Photon spectra of all four models integrated over the 
region from the CD to the FS. Upper three panels: Models B to 
D are compared to Model A and are split into individual compo- 
nents for different emission mechanisms: 7r°-decay (red), IC (blue), 
bremsstrahlung (green) and synchrotron radiation (black). Solid 
lines represent spectra for Model B, C and D in each panel while 
Model A is shown as dashed lines. Bottom panel: The contribu- 
tions from all mechanisms are summed for each model: Model A 
(thick solid), Model B (thin solid), Model C (dash) and Model D 
(dash-dot). 

cient Z?Koi °c P 1 ^ 3 , hardening the spectrum even more, 
as shown in Fig. |4j The photon spectra reflect the hard 
particle spectra. 

With a number density of timc = 10 3 cm~ 3 , and a 
column density of ncoi ~ 10 21 " 22 cm~ 2 in the 'molecular 
cloud,' 7r°-decay is the main 7-ray source for all models, 
followed by relativistic bremsstrahlung and then IC emis- 
sion, as shown in Fig. [SI For Models A, B, C and D, the 
separation between the FS and the inner edge of the MC 



Fig. 6. — Same as Fig. [5] but the photon spectra are now inte- 
grated over the shell molecular cloud volume. The strong momen- 
tum dependent diffusion results in hard spectra and a change in the 
relative dominance between the emission mechanisms, compared to 
spectra integrated over the SNR. 



is found to be around 4.3, 3.8, 2.3 and 5.9 pc respectively 
at t = 1000 yr. 

For the environmental parameters studied here, the 
emission at all wavebands from the molecular cloud is 
weaker than that from the SNR shell, but the difference 
depends on the photon energy. With the assumption 
Bmc = 3 /iG, the X-ray synchrotron flux from the molec- 
ular cloud stays at the ISM level and will be difficult to 
detect. The GeV 7-ray flux is more model-dependent and 
the flux stays around a factor of 10-100 smaller than the 
flux from the SNR. For the TcV flux, which is detected by 
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Fig. 7. — Photon flux maps projected along the line-of-sight for Model A (left panels) and Model B (right panels), with a hemisphere 
molecular cloud centered at pixel coordinate (50,25,25) with a radius of 3.2 pc. The horizontal and virtical scales are in pixels where the 
pixel size is 0.38 pc X 0.38 pc. The upper panels are integrated over the energy range 1 — 300 GeV, while the bottom panels are integrated 
gies £ 7 > 1 TeV. The color scale is logarithmic in log(A r 7 /cm 2 /s). The dashed circle in each panel indicates the position of the 



over energies 



contact discontinuity. 

atmospheric Cherenkov telescopes, the molecular cloud 
emission is about 1 — 10% of that from the SNR. 



4.4. Boardband Images and Projected Emission Profiles 

Multi-wavelength projection maps are useful for 
studying the energy-dependent morphology of SNRs. We 
use our hemispherical 'molecular cloud' example (Fig.[T]) 
to calculate 2-D projection maps in various energy bands 
at tsNR = 1000 yr for Models A and B. After the pho- 
ton emissivity is calculated in each boxel in the 3-D 
simulation box, we perform a line-of-sight projection 
through the box. We choose four energy bands: (i) 
soft X-rays with E 1 = 1 — 5 keV; (ii) hard X-rays with 
E 1 = 5 - lOkcV; (hi) £7 7 = 1 - 300 GeV; and (iv) 
E 1 > 1 TeV. The parameters we use result in a col- 
umn density of ~ 10 21 cm -2 for the cloud, which is small 
enough to ignore in the present context. 

Fig. [7] shows the 7-ray projected flux maps in 



log(A 7 /cm 2 /s) at a source distance of -Dsnr = lkpc for 
Models A (NL) and B (TP) in the GeV and TeV bands 
[i.e., bands (hi) and (iv)]. The color scales are different 
for the GeV and TeV images, but the spatial resolution 
is the same. The difference between the example with ef- 
ficient DSA (left panels) and the test-particle case (right 
panels) is mainly one of intensity if only GeV- TeV emis- 
sion is concerned. In both cases, the brightest regions 
of the SNR are considerably brighter than the cloud and 
for the test-particle case (Model B), the cloud is almost 
invisible on these scales. For the SNR in both the NL 
and TP cases, the region between the CD and FS clearly 
shows up in the maps even with the projection through 
the remnant. 8 There is also a clear limb darkening effect 
from projection seen at the edge of both remnants and at 
the edge of the cloud in the NL case. At the molecular 

8 The dashed circle in each panel shows the position of the CD 
at 1000 yr. 
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Fig. 8. — Line-of-sight emission profiles in the radial direction 
from the SNR center at R = to the hemispheric 'molecular cloud' 
centered at R = 10 pc. Results for Model A (solid curves) and 
Model B (dashed curves) are displayed. Markers represent the 
centers of the cubic spatial bins which are 0.385 pc on a side. Four 
wavebands are considered - (i) Soft X-rays from 1 to 5keV (circle); 
(ii) Hard X-rays from 5 to lOkeV (triangle); (iii) 1 to 300 GeV 
7-rays (square) and (iv) 7-rays with energy above 1 TeV (cross) . 
The vertical axis shows the photon flux in log(A r 7 /cm 2 /s) for each 
pixel bin. 



cloud, however, there is no noticable drop in intensity 
towards the center of the cloud as occurs for the SNR. 

These details, of course, depend on the particular 
parameters we have chosen but some general statements 
can be made. Unless there is a source of soft photons 
associated with the external material, the brightness of 
the external material (MC) compared to the SNR, the 
Imc/Isnr ratio, will be independent of the density ratio 
n-MC /n p if IC dominates the GeV- TeV emission. If pion- 
decay or bremsstrahlung dominate, the /mc/^snr ratio 
will scale approximately as the first power of the den- 
sity ratio, riMc/ n p- I n cases, /mc/^SNr will decrease 
with the distance the external material is from the FS. 
Another important result, which is implicit in Fig. [7]and 
important for comparing pion-decay and IC emission, is 
that emission from the SNR and the external material 
must be considered together. To first order, an increase 
in acceleration efficiency or ambient matter density not 
solely associated with the cloud, n p , will leave /mc/^snr 
unchanged. 

In Fig. [8] we show emission calculated along a hori- 
zontal line from the center of the remnant at R = across 
the molecular cloud for all four energy bands. These 
fluxes are determined, as are the 2-D maps, by summing 
the emission from each boxel along a line-of-sight. The 
plateaus on the left hand side of the plot within R < 5 pc 
show emission from the SNR. The subtle increase of the 
projected flux with R in this region is the result of pro- 
jection through the shell of material between the CD and 
the FS. Beyond R ~ 5 pc, the fluxes drop abruptly to the 
ISM level. Here, escaping CRs stream through the ISM 
with a large diffusion coefficient DkoI- At R ~ 6.8 pc, 
the CRs impact the hemisphere 'molecular cloud' with 
riMC = 10 3 cm -3 and the fluxes for energy bins (iii) and 
(iv) increase by almost 2 orders of magnitude from the 
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Fig. 9. — Same as Fig. [8] with individual emission processes 
shown, i.e., pion-decay (circles), IC (triangles), and bremsstrahlung 
(squares). Synchrotron emission has a negligible flux in the 7-ray 
energy band. 



ISM level. These photons are mainly from pion-decay. 
There is no increase at the edge of the cloud for energy 
bins (i) and (ii) since this emission is totally from syn- 
chrotron and we have assumed the field in the molecular 
cloud equals the ISM field, Bmc = -Bism- 

Fig. [5] shows the emission profiles in bands (iii) and 
(iv) for Models A and B separated into individual emis- 
sion mechanisms. While the total fluxes at these energies 
depend strongly on acceleration efficiency, the Tmc / ^snr 
ratio is much less sensitive to e acc , as mentioned above. 
There is no increase in IC emission at the edge of the 
cloud near R ~ 6.9 pc since we only consider electron 
scattering off the CMB. Unless there is an additional 
source of photons associated with the external material, 
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IC will be strongly suppressed relative to pion-decay in 
external material. 

5. SUMMARY AND CONCLUSIONS 

We have introduced a 3-D simulation of an evolv- 
ing SNR where the nonlinear acceleration of CRs is cou- 
pled to the SNR evolution. The model follows the diffu- 
sion and interaction of CRs within the spherically sym- 
metric remnant, as well as high-energy CRs that escape 
from the forward shock and diffuse into the surrounding 
medium. For any set of model and environmental pa- 
rameters, and for arbitrary distributions of matter sur- 
rounding the remnant, wc can calculate broadband pho- 
ton spectra and obtain line-of-sight projections and mor- 
phologies that will allow for efficient comparisons with 
observations in various energy bands. 

We have illustrated the capabilities of this simula- 
tion with several models that differ from each other in 
the CR acceleration efficiency, e acc , the ambient ISM pro- 
ton density, n p , and the matter distribution of a 'molec- 
ular cloud' external to the SNR. Of course, all of the re- 
sults discussed here assume particular values for param- 
eters, such as a shocked electron to proton temperature 
ratio T e /T p = 1 and an electron/proton ratio at rela- 
tivistic energies (e/p) Ie \ = 0.01. These parameters are 
critical for understanding DSA and applying the mech- 
anism to astrophysical sources yet they are poorly con- 
strained by both observations and theory. For instance, 
the value of T e /T p determines the importance of brcms- 
strahlung compared to synchrotron in the X-ray range 
and also strongly influe nces the thermal X-ray line spec- 
tra (|Ellison et al.ll2007f) . The (e/p) rc i ratio is the most 
important factor after the ambient density determining 
the relative intensity of IC and pion-decay emission at 
GeV-TeV energies. The confirmation of CR ion produc- 
tion in SNRs depends on this parameter. 

Other important parameters of DSA that remain 
uncertain are the injection and acceleration efficiencies, 
the amount of magnetic compression and amplification 
that occurs, and the diffusion coefficient of escaping par- 
ticles as they leave the shock, which must differ sub- 
stantially from Bohm diffu sion fe.g.. iBlasi et al.1 l2007t 
lEllison fe Vladimirovll2"008t ). Due to the still limited dy- 
namic range of particle-in-ccll simulations, and the lack 
of strong, nonlinear shocks producing relativistic par- 
ticles in the hcliosphcre, we believe young SNRs are 
the best 'laboratory' for studying NL-DSA. Broadband 
observations matched against self-consistent nonlinear 
models currently provide the best constraints on these 
important parameters. 

There are three important aspects of our 3-D simu- 
lation that are new and extend the large body of existing 
work on DSA in SNRs. One is that the simulation con- 
sistently models high-energy CRs that escape from the 
forward shock of the SNR with the evolution of the SNR 



itself. In NL-DSA, the fraction of total explosion energy 
that ends up in escaping particles can be large (see Ta- 
ble [1]) and we believe this is the first work to include 
these particles in a coherent emission model. Second, 
the 3-D simulation box allows for the modeling of CR in- 
teractions in arbitrary mass distributions outside of the 
SNR. This feature is essential for producing 2-D projec- 
tion maps that can be compared with current and future 
observations. These maps, tuned to match the instru- 
ment response of telescopes, will serve to help determine 
the importance of pre-SN shells and/or nearby molecular 
clouds in producing 7-ray emission. Third, the simula- 
tion platform is extremely flexible making it straightfor- 
ward to add important effects not present in this pre- 
liminary model. These generalizations include shock ac- 
celeration and heating at the reverse shock as well as 
the forward shock, pre-SN winds for Type II SNe, var- 
ious forms for particle diffusion in the ISM, production 
and interaction of heavy CR ions, and a parameterized 
representation of magnetic field amplification. 

Another physical effect that may importantly 
influence the photon spectrum is anisotropy from 
angular-dependent interactions. These include an 
angular-dependent neutral pion production cross-section 
(jKarlsson fc Kamad |2007( ) and anisotropic IC scat- 
tering with photon field s other than the CMB 
(|Moskalenko fc Strong! [20001 . Preliminary results show 
that anisotropies can change the spectral shape and flux 
of the observed photons drastically. When anisotropic 
interactions are implemented, the projection maps we 
calculate will show how the observed flux depends on 
the orientation of the FS and molecular cloud with re- 
spect to the line-of-sight. We leave this issue to future 
studies. 

Finally, in addition to modeling the photon emis- 
sion from SNRs, our model can also determine the total 
contribution of CR ions and electrons injected into the 
Galaxy from an individual SNe over its lifetime. This can 
serve as input to Ga laxy-scale propaga tion models (for 
example, GALPROP. IStrong etaDl2007h and also add to 
our knowledge on the Galactic 7-ray diffuse emission. 
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